home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / b / b.lha / B / src / bint / b1num.c < prev    next >
Encoding:
C/C++ Source or Header  |  1988-11-24  |  6.9 KB  |  314 lines
  1. /* Copyright (c) Stichting Mathematisch Centrum, Amsterdam, 1985. */
  2.  
  3. /*
  4.   $Header: b1num.c,v 1.4 85/08/22 16:51:59 timo Exp $
  5. */
  6.  
  7. /* B numbers, basic external interface */
  8.  
  9. #include "b.h"
  10. #include "b0con.h"
  11. #include "b1obj.h"
  12. #include "b1num.h"
  13.  
  14. /*
  15.  * This file contains operations on numbers that are not predefined
  16.  * B functions (the latter are defined in `bfun.c').
  17.  * This includes conversion of numeric B values to C `int's and
  18.  * `double's, (numval() and intval()),
  19.  * but also utilities for comparing numeric values and hashing a numeric
  20.  * value to something usable as initialization for the random generator
  21.  * without chance of overflow (so numval(v) is unusable).
  22.  * It is also possible to build numbers of all types using mk_integer,
  23.  * mk_exact (or mk_rat) and mk_approx.  Note the rather irregular
  24.  * (historical!) argument structure for these: mk_approx has a
  25.  * `double' argument, where mk_exact and mk_rat have two values
  26.  * which must be `integer' (not `int's!) as arguments.
  27.  * The random number generator used by the DRAW and CHOOSE statements
  28.  * is also in this file.
  29.  */
  30.  
  31. /*
  32.  * ival is used internally by intval() and large().
  33.  * It converts an integer to double precision, yielding
  34.  * a huge value when overflow occurs (but giving no error).
  35.  */
  36.  
  37. Hidden double ival(u) integer u; {
  38.     double x = 0;
  39.     register int i;
  40.  
  41.     if (IsSmallInt(u)) return SmallIntVal(u);
  42.     for (i = Length(u)-1; i >= 0; --i) {
  43.         if (x >= Maxreal/BASE)
  44.             return Msd(u) < 0 ? -Maxreal : Maxreal;
  45.         x = x*BASE + Digit(u, i);
  46.     }
  47.  
  48.     return x;
  49. }
  50.  
  51.  
  52. /* Determine if a value may be converted to an int */
  53.  
  54. Visible bool large(v) value v; {
  55.     double r;
  56.     if (!Is_number(v) || !integral(v)) {
  57.         error(MESS(1300, "number not an integer"));
  58.         return No;
  59.     }
  60.     if (Rational(v)) v= (value) Numerator((rational)v);
  61.     r= ival((integer)v);
  62.     if (r > Maxint) return Yes;
  63.     return No;
  64. }
  65.  
  66. /* return the C `int' value of a B numeric value, if it exists. */
  67.  
  68. Visible int intval(v) value v; {
  69.     /* v must be an Integral number or a Rational with Denominator==1
  70.         which may result from n round x [via mk_exact]!. */
  71.     double i;
  72.     if (IsSmallInt(v)) i= SmallIntVal(v);
  73.     else {
  74.         if (!Is_number(v)) syserr(MESS(1301, "intval on non-number"));
  75.         if (!integral(v)) {
  76.             error(MESS(1302, "number not an integer"));
  77.             return 0;
  78.         }
  79.         if (Rational(v)) v= (value) Numerator((rational)v);
  80.         i= ival((integer)v);
  81.     }
  82.     if (i > Maxint || i < -Maxint) {
  83.         error(MESS(1303, "exceedingly large integer"));
  84.         return 0;
  85.     }
  86.     return (int) i;
  87. }
  88.  
  89.  
  90. /* convert an int to an intlet */
  91.  
  92. Visible intlet propintlet(i) int i; {
  93.     if (i > Maxintlet || i < -Maxintlet) {
  94.         error(MESS(1304, "exceedingly large integer"));
  95.         return 0;
  96.     }
  97.     return i;
  98. }
  99.  
  100.  
  101. /*
  102.  * determine if a number is integer
  103.  */
  104.  
  105. Visible bool integral(v) value v; {
  106.     if (Integral(v) || (Rational(v) && Denominator((rational)v) == int_1))
  107.         return Yes;
  108.     else return No;
  109. }
  110.  
  111. /*
  112.  * mk_exact makes an exact number out of two integers.
  113.  * The third argument is the desired number of digits after the decimal
  114.  * point when the number is printed (see comments in `bfun.c' for
  115.  * `round' and code in `bnuC.c').
  116.  * This printing size (if positive) is hidden in the otherwise nearly
  117.  * unused * 'Size' field of the value struct
  118.  * (cf. definition of Rational(v) etc.).
  119.  */
  120.  
  121. Visible value mk_exact(p, q, len) integer p, q; int len; {
  122.     rational r = mk_rat(p, q, len);
  123.  
  124.     if (Denominator(r) == int_1 && len <= 0) {
  125.         integer i = (integer) Copy(Numerator(r));
  126.         release((value) r);
  127.         return (value) i;
  128.     }
  129.  
  130.     return (value) r;
  131. }
  132.  
  133. #define uSMALL 1
  134. #define uINT 2
  135. #define uRAT 4
  136. #define uAPP 8
  137. #define vSMALL 16
  138. #define vINT 32
  139. #define vRAT 64
  140. #define vAPP 128
  141.  
  142. Visible relation numcomp(u, v) value u, v; {
  143.     int tu, tv; relation s;
  144.  
  145.     if (IsSmallInt(u)) tu = uSMALL;
  146.     else if (Length(u) >= 0) tu = uINT;
  147.     else if (Length(u) <= -2) tu = uRAT;
  148.     else tu = uAPP;
  149.     if (IsSmallInt(v)) tv = vSMALL;
  150.     else if (Length(v) >= 0) tv = vINT;
  151.     else if (Length(v) <= -2) tv = vRAT;
  152.     else tv = vAPP;
  153.  
  154.     switch (tu|tv) {
  155.  
  156.     case uSMALL|vSMALL: return SmallIntVal(u) - SmallIntVal(v);
  157.  
  158.     case uSMALL|vINT:
  159.     case uINT|vSMALL:
  160.     case uINT|vINT: return int_comp((integer)u, (integer)v);
  161.  
  162.     case uSMALL|vRAT:
  163.     case uINT|vRAT:
  164.         u = (value) mk_rat((integer)u, int_1, 0);
  165.         s = rat_comp((rational)u, (rational)v);
  166.         release(u);
  167.         return s;
  168.  
  169.     case uRAT|vRAT:
  170.         return rat_comp((rational)u, (rational)v);
  171.  
  172.     case uRAT|vSMALL:
  173.     case uRAT|vINT:
  174.         v = (value) mk_rat((integer)v, int_1, 0);
  175.         s = rat_comp((rational)u, (rational)v);
  176.         release(v);
  177.         return s;
  178.  
  179.     case uSMALL|vAPP:
  180.     case uINT|vAPP:
  181.     case uRAT|vAPP:
  182.         u = approximate(u);
  183.         s = app_comp((real)u, (real)v);
  184.         release(u);
  185.         return s == 0 ? -1 : s; /* u < ~u = v */
  186.  
  187.     case uAPP|vAPP:
  188.         return app_comp((real)u, (real)v);
  189.  
  190.     case uAPP|vSMALL:
  191.     case uAPP|vINT:
  192.     case uAPP|vRAT:
  193.         v = approximate(v);
  194.         s = app_comp((real)u, (real)v);
  195.         release(v);
  196.         return s == 0 ? 1 : s; /* u = ~v > v */
  197.  
  198.     default: syserr(MESS(1305, "num_comp")); /* NOTREACHED */
  199.  
  200.     }
  201. }
  202.  
  203.  
  204. /*
  205.  * Deliver 10**n, where n is a (maybe negative!) C integer.
  206.  * The result is a value (integer or rational, actually).
  207.  */
  208.  
  209. Visible value tento(n) int n; {
  210.     if (n < 0) {
  211.         integer i= int_tento(-n);
  212.         value v= (value) mk_exact(int_1, i, 0);
  213.         release((value) i);
  214.         return v;
  215.     }
  216.     return (value) int_tento(n);
  217. }
  218.  
  219.  
  220. /*
  221.  * numval returns the numerical value of any numeric B value
  222.  * as a C `double'.
  223.  */
  224.  
  225. Visible double numval(u) value u; {
  226.     double expo, frac;
  227.  
  228.     if (!Is_number(u)) {
  229.         error(MESS(1306, "value not a number"));
  230.         return 0.0;
  231.     }
  232.     u = approximate(u);
  233.     expo = Expo((real)u), frac = Frac((real)u);
  234.     release(u);
  235.     if (expo > Maxexpo) {
  236.         error(MESS(1307, "approximate number too large to be handled"));
  237.         return 0.0;
  238.     }
  239.     if(expo < Minexpo) return 0.0;
  240.     return ldexp(frac, (int)expo);
  241. }
  242.  
  243.  
  244. /*
  245.  * Random numbers
  246.  */
  247.  
  248.  
  249. /*
  250.  * numhash produces a `double' number that depends on the value of
  251.  * v, useful for initializing the random generator.
  252.  * Needs rewriting, so that for instance numhash(n) never equals n.
  253.  * The magic numbers here are chosen at random.
  254.  */
  255.  
  256. Visible double numhash(v) value v; {
  257.     if (Integral(v)) {
  258.         double d = 0;
  259.         register int i;
  260.  
  261.         if (IsSmallInt(v)) return SmallIntVal(v);
  262.         for (i = Length(v) - 1; i >= 0; --i) {
  263.             d *= 2;
  264.             d += Digit((integer)v, i);
  265.         }
  266.  
  267.         return d;
  268.     }
  269.  
  270.     if (Rational(v))
  271.         return .777 * numhash((value) Numerator((rational)v)) +
  272.                .123 * numhash((value) Denominator((rational)v));
  273.  
  274.     return numval(v); /* Which fails for HUGE reals.  Never mind. */
  275. }
  276.  
  277.  
  278. /* Initialize the random generator */
  279.  
  280. double lastran;
  281.  
  282. Hidden Procedure setran (seed) double seed; {
  283.     double x;
  284.  
  285.     x = seed >= 0 ? seed : -seed;
  286.     while (x >= 1) x /= 10;
  287.     lastran = floor(67108864.0*x);
  288. }
  289.  
  290. Visible Procedure set_random(v) value v; {
  291.     setran((double) hash(v));
  292. }
  293.  
  294.  
  295. /* Return a random number in [0, 1). */
  296.  
  297. Visible value random() {
  298.     double p;
  299.  
  300.     p = 26353589.0 * lastran + 1;
  301.     lastran = p - 67108864.0*floor(p/67108864.0);
  302.  
  303.     return (value) mk_approx(lastran / 67108864.0, 0.0);
  304. }
  305.  
  306. Visible Procedure initnum() {
  307.     rat_init();
  308.     setran((double) SEED);
  309. }
  310.  
  311. Visible Procedure endnum() {
  312.     endrat();
  313. }
  314.